Differential Expression Analysis
Copyright © 2024 Daniel Zimmermann. daniel.zimmermann.dz@outlook.com
You may use, distribute and modify this code under the terms of the MIT license.
You should have received a copy of the MIT license with this file. If not, please visit: https://opensource.org/license/mit
Setup
source("helpers.R")## Lade nötiges Paket: S4Vectors## Lade nötiges Paket: stats4## Lade nötiges Paket: BiocGenerics##
## Attache Paket: 'BiocGenerics'## Die folgenden Objekte sind maskiert von 'package:stats':
##
## IQR, mad, sd, var, xtabs## Die folgenden Objekte sind maskiert von 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
## table, tapply, union, unique, unsplit, which.max, which.min##
## Attache Paket: 'S4Vectors'## Das folgende Objekt ist maskiert 'package:utils':
##
## findMatches## Die folgenden Objekte sind maskiert von 'package:base':
##
## expand.grid, I, unname## Lade nötiges Paket: IRanges##
## Attache Paket: 'IRanges'## Das folgende Objekt ist maskiert 'package:grDevices':
##
## windows## Lade nötiges Paket: GenomicRanges## Lade nötiges Paket: GenomeInfoDb## Lade nötiges Paket: SummarizedExperiment## Lade nötiges Paket: MatrixGenerics## Lade nötiges Paket: matrixStats##
## Attache Paket: 'MatrixGenerics'## Die folgenden Objekte sind maskiert von 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Lade nötiges Paket: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attache Paket: 'Biobase'## Das folgende Objekt ist maskiert 'package:MatrixGenerics':
##
## rowMedians## Die folgenden Objekte sind maskiert von 'package:matrixStats':
##
## anyMissing, rowMedians##
## Attache Paket: 'dplyr'## Das folgende Objekt ist maskiert 'package:Biobase':
##
## combine## Das folgende Objekt ist maskiert 'package:matrixStats':
##
## count## Die folgenden Objekte sind maskiert von 'package:GenomicRanges':
##
## intersect, setdiff, union## Das folgende Objekt ist maskiert 'package:GenomeInfoDb':
##
## intersect## Die folgenden Objekte sind maskiert von 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union## Die folgenden Objekte sind maskiert von 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union## Die folgenden Objekte sind maskiert von 'package:BiocGenerics':
##
## combine, intersect, setdiff, union## Die folgenden Objekte sind maskiert von 'package:stats':
##
## filter, lag## Die folgenden Objekte sind maskiert von 'package:base':
##
## intersect, setdiff, setequal, union##
## Attache Paket: 'magrittr'## Das folgende Objekt ist maskiert 'package:GenomicRanges':
##
## subtract## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%within%() masks IRanges::%within%()
## ✖ dplyr::collapse() masks IRanges::collapse()
## ✖ dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count() masks matrixStats::count()
## ✖ dplyr::desc() masks IRanges::desc()
## ✖ tidyr::expand() masks S4Vectors::expand()
## ✖ tidyr::extract() masks magrittr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks S4Vectors::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename() masks S4Vectors::rename()
## ✖ lubridate::second() masks S4Vectors::second()
## ✖ lubridate::second<-() masks S4Vectors::second<-()
## ✖ purrr::set_names() masks magrittr::set_names()
## ✖ dplyr::slice() masks IRanges::slice()
## ✖ magrittr::subtract() masks GenomicRanges::subtract()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Lade nötiges Paket: AnnotationDbi
##
##
## Attache Paket: 'AnnotationDbi'
##
##
## Das folgende Objekt ist maskiert 'package:dplyr':
##
## select
##
##
##
##
## Lade nötiges Paket: data.table
##
##
## Attache Paket: 'data.table'
##
##
## Die folgenden Objekte sind maskiert von 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
##
##
## Das folgende Objekt ist maskiert 'package:purrr':
##
## transpose
##
##
## Die folgenden Objekte sind maskiert von 'package:dplyr':
##
## between, first, last
##
##
## Das folgende Objekt ist maskiert 'package:SummarizedExperiment':
##
## shift
##
##
## Das folgende Objekt ist maskiert 'package:GenomicRanges':
##
## shift
##
##
## Das folgende Objekt ist maskiert 'package:IRanges':
##
## shift
##
##
## Die folgenden Objekte sind maskiert von 'package:S4Vectors':
##
## first, second
##
##
## Lade nötiges Paket: grid
##
## ========================================
## ComplexHeatmap version 2.22.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
##
##
## ========================================
## circlize version 0.4.16
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================
##
##
##
## Attache Paket: 'plotly'
##
##
## Das folgende Objekt ist maskiert 'package:ComplexHeatmap':
##
## add_heatmap
##
##
## Das folgende Objekt ist maskiert 'package:AnnotationDbi':
##
## select
##
##
## Das folgende Objekt ist maskiert 'package:ggplot2':
##
## last_plot
##
##
## Das folgende Objekt ist maskiert 'package:IRanges':
##
## slice
##
##
## Das folgende Objekt ist maskiert 'package:S4Vectors':
##
## rename
##
##
## Das folgende Objekt ist maskiert 'package:stats':
##
## filter
##
##
## Das folgende Objekt ist maskiert 'package:graphics':
##
## layoutdir.create("results/DESeq", showWarnings = F, recursive = T)
knitr::opts_chunk$set(fig.width = 10, dpi = 300, results = "hold", fig.show = "hold")## Counts plot
countsToPlot <- c("CD274", "BRAF", "BCL2", "VEGFC", "PDGFC")
## Heatmaps
hm_scale_by_row <- TRUE
hm_max_rows <- 50
hmap_colors <- colorRamp2(c(-1.5, 0, 1.5), c("darkblue", "white", "darkred"))
annotation_colors <- c("nHEM" = "#F8766D",
"MCM1G" = "#7CAE00",
"MCMDLN" = "#00BFC4",
"MCMDLN_3D" = "#C77CFF")
## Volcano plot
vp_max_labels <- 0
vp_lfc_limit <- NAImport data
# Read metadata
samples <- read.delim("data/metadata.txt") %>% as_tibble()
# Read featurecounts output
countdata <- read_counts(samples$Files, "data/featurecounts/", sample_IDs = samples$Sample_ID)
# Read GTF file
genedata <- read_genes("data/genes.gtf", countdata$Geneid)## Automatically detected number of rows to skip: 0
## List of features in column 9:
## -----------------------------
## gene_biotype
## gene_id
## gene_name
## gene_source
## transcript_id
## transcript_name
## transcript_source
## tss_idPrepare data for analysis
countdata.filtered <- filter_counts(countdata)
# Move gene info to separate dataframe and remove from count data
genedata.filtered <- genedata %>%
filter(ENSEMBL %in% rownames(countdata.filtered))
samples %<>%
mutate(Type = str_remove(Sample_ID, "(_[23]D)?_\\d$"),
Class = factor(Class, c("nHEM", "MCM1G", "MCMDLN", "MCMDLN_3D")))
summary(samples)## Sample_ID Files Replicate Annotation
## Length:12 Length:12 Min. :1 Length:12
## Class :character Class :character 1st Qu.:1 Class :character
## Mode :character Mode :character Median :2 Mode :character
## Mean :2
## 3rd Qu.:3
## Max. :3
## Type Dimension Class
## Length:12 Length:12 nHEM :3
## Class :character Class :character MCM1G :3
## Mode :character Mode :character MCMDLN :3
## MCMDLN_3D:3
##
## Run DESeq analysis
# Create DESeq dataset
dds <- DESeqDataSetFromMatrix(countData = countdata.filtered, colData = samples, design = ~Class)
# Run DESeq analysis
dds = DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingnormalized_counts = counts(dds, normalized = T)Overview and QC plots
# Size Factors for Normalization
sizeFactors(dds)
# Dispersion plot for QC
plotDispEsts(dds)
# transform data with regularized log for visualization
vsd <- vst(dds, blind = T)
# Principal component analysis
plot_pca(vsd)
# Correlation heatmap
samplecorr <- cor(assay(vsd))
rownames(samplecorr) <- samples$Sample_ID
hmap_annotation <- samples %>%
dplyr::select(c(Class))
Heatmap(samplecorr, c("darkblue", "khaki1", "darkred"), name = "Correlation",
show_column_names = FALSE, row_dend_reorder = FALSE, column_dend_reorder = FALSE,
top_annotation = HeatmapAnnotation(Class = hmap_annotation$Class,
col = list(Class = annotation_colors),
show_annotation_name = FALSE))## nHEM_1 nHEM_2 nHEM_3 MCM1G_1 MCM1G_2 MCM1G_3
## 0.4848247 0.6124169 0.7109039 1.4445276 1.3167976 1.1891653
## MCMDLN_1 MCMDLN_2 MCMDLN_3 MCMDLN_3D_1 MCMDLN_3D_2 MCMDLN_3D_3
## 0.9962995 1.0757996 1.2974828 1.1285081 1.0287994 1.5023460
Extract results from DESeq object
res.1GvsHEM = results(dds, contrast = c("Class", "MCM1G", "nHEM"), alpha = 0.05)
res.DLNvsHEM = results(dds, contrast = c("Class", "MCMDLN", "nHEM"), alpha = 0.05)
res.DLNvs1G = results(dds, contrast = c("Class", "MCMDLN", "MCM1G"), alpha = 0.05)
res.3Dvs2D = results(dds, contrast = c("Class", "MCMDLN_3D", "MCMDLN"), alpha = 0.05)
plotMA(res.1GvsHEM, ylim=c(-2, 2), main = "MCM1G vs nHEM")
plotMA(res.DLNvsHEM, ylim=c(-2, 2), main = "MCMDLN vs nHEM")
plotMA(res.DLNvs1G, ylim=c(-2, 2), main = "MCMDLN vs MCM1G")
plotMA(res.3Dvs2D, ylim=c(-2, 2), main = "MCMDLN (3D) vs MCMDLN (2D)")
print("MCM1G vs nHEM")
summary(res.1GvsHEM)
print("MCMDLN vs nHEM")
summary(res.DLNvsHEM)
print("MCMDLN vs MCM1G")
summary(res.DLNvs1G)
print("MCMDLN (3D) vs MCMDLN (2D)")
summary(res.3Dvs2D)## [1] "MCM1G vs nHEM"
##
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 9536, 21%
## LFC < 0 (down) : 11374, 25%
## outliers [1] : 13, 0.028%
## low counts [2] : 9896, 21%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "MCMDLN vs nHEM"
##
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 10379, 22%
## LFC < 0 (down) : 11900, 26%
## outliers [1] : 13, 0.028%
## low counts [2] : 9896, 21%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "MCMDLN vs MCM1G"
##
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 8332, 18%
## LFC < 0 (down) : 7563, 16%
## outliers [1] : 13, 0.028%
## low counts [2] : 10795, 23%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "MCMDLN (3D) vs MCMDLN (2D)"
##
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 5219, 11%
## LFC < 0 (down) : 5392, 12%
## outliers [1] : 13, 0.028%
## low counts [2] : 15293, 33%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
Perform shrinkage
res.1GvsHEM <- lfcShrink(dds = dds, res = res.1GvsHEM, type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041res.DLNvsHEM <- lfcShrink(dds = dds, res = res.DLNvsHEM, type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041res.DLNvs1G <- lfcShrink(dds = dds, res = res.DLNvs1G, type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041res.3Dvs2D <- lfcShrink(dds = dds, res = res.3Dvs2D, type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041plotMA(res.1GvsHEM, ylim=c(-2, 2), main = "MCM1G vs nHEM")
plotMA(res.DLNvsHEM, ylim=c(-2, 2), main = "MCMDLN vs nHEM")
plotMA(res.DLNvs1G, ylim=c(-2, 2), main = "MCMDLN vs MCM1G")
plotMA(res.3Dvs2D, ylim=c(-2, 2), main = "MCMDLN (3D) vs MCMDLN (2D)")
Results
results_all <- list(MCM1G_vs_nHEM = res.1GvsHEM,
MCMDLN_vs_nHEM = res.DLNvsHEM,
MCMDLN_vs_MCM1G = res.DLNvs1G,
MCMDLN_3D_vs_2D = res.3Dvs2D)
results_all %>%
lapply(function(x) {
x %>%
as.data.frame() %>%
filter(padj < 0.05 & abs(log2FoldChange) >= 1) %>%
dplyr::select(log2FoldChange) %>%
sign() %>%
mutate(Expression = ifelse(log2FoldChange == 1, "Overexpressed", "Underexpressed")) %>%
dplyr::select(!log2FoldChange) %>%
group_by(Expression) %>%
count()
}) %>%
bind_rows(.id = "Comparison") %>%
pivot_wider(names_from = Expression, values_from = n) %>%
mutate(Total = Overexpressed + Underexpressed)Expression levels (all samples)
counts.plot <- lapply(countsToPlot, function(x) {
plotCounts(dds, genedata.filtered[genedata.filtered$Name == x,]$ENSEMBL,
intgroup = c("Class"), returnData = TRUE)
}) %>%
set_names(countsToPlot) %>%
lapply(rownames_to_column, var = "Sample_ID") %>%
bind_rows(.id = "Gene") %>%
group_by(Gene, Class) %>%
summarise(mean_count = mean(count), sd_count = sd(count), .groups = "drop") %>%
mutate(min_error = ifelse(mean_count - sd_count >= 0, mean_count - sd_count, 0),
max_error = mean_count + sd_count)
ggplot(counts.plot, aes(x = Gene, y = mean_count, fill = Class)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = min_error, ymax = max_error),
position = position_dodge(0.9), width = 0.2) +
labs(x = "", y = "Mean Count", fill = "") +
theme_pubr() +
theme(legend.position = "right", legend.text = element_text(size = 12), panel.grid.major.y = element_line(linewidth = 0.5))
MCM1G vs nHEM
res.1GvsHEM.tb <- res.1GvsHEM %>%
data.frame() %>%
rownames_to_column(var = "ENSEMBL") %>%
as_tibble() %>%
left_join(genedata.filtered, ., by = "ENSEMBL") %>%
mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
res.1GvsHEM.sig <- res.1GvsHEM.tb %>%
filter(threshold == TRUE)
pvalues <- res.1GvsHEM.sig$padj
names(pvalues) <- rownames(res.1GvsHEM.sig)
norm.1GvsHEM.sig <- normalized_counts %>%
data.frame() %>%
dplyr::select(samples$Sample_ID[samples$Class %in% c("MCM1G", "nHEM")]) %>%
filter(rownames(.) %in% res.1GvsHEM.sig$ENSEMBL)
write.csv(res.1GvsHEM.tb %>% dplyr::select(-threshold),
file = "results/DESeq/MCM1G_vs_nHEM_unfiltered.csv",
row.names = FALSE)
write.csv(res.1GvsHEM.sig %>% dplyr::select(-threshold),
file = "results/DESeq/MCM1G_vs_nHEM_filtered.csv",
row.names = FALSE)plot_heatmap(norm.1GvsHEM.sig, hmap_colors, annotation_colors,
genedata.filtered, samples, pvalues, max_rows = hm_max_rows)
volcanoplot(res.1GvsHEM.tb, title = "MCM1G vs. nHEM", max_labels = vp_max_labels,
lfc_limit = vp_lfc_limit)
volcanoplot_interactive(res.1GvsHEM.tb, title = "MCM1G vs. nHEM")MCMDLN vs nHEM
res.DLNvsHEM.tb <- res.DLNvsHEM %>%
data.frame() %>%
rownames_to_column(var = "ENSEMBL") %>%
as_tibble() %>%
left_join(genedata.filtered, ., by = "ENSEMBL") %>%
mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
res.DLNvsHEM.sig <- res.DLNvsHEM.tb %>%
filter(threshold == TRUE)
pvalues <- res.DLNvsHEM.sig$padj
names(pvalues) <- rownames(res.DLNvsHEM.sig)
norm.DLNvsHEM.sig <- normalized_counts %>%
data.frame() %>%
dplyr::select(samples$Sample_ID[samples$Class %in% c("MCMDLN", "nHEM")]) %>%
filter(rownames(.) %in% res.DLNvsHEM.sig$ENSEMBL)
write.csv(res.DLNvsHEM.tb %>% dplyr::select(-threshold),
file = "results/DESeq/MCMDLN_vs_nHEM_unfiltered.csv",
row.names = FALSE)
write.csv(res.DLNvsHEM.sig %>% dplyr::select(-threshold),
file = "results/DESeq/MCMDLN_vs_nHEM_filtered.csv",
row.names = FALSE)plot_heatmap(norm.DLNvsHEM.sig, hmap_colors, annotation_colors,
genedata.filtered, samples, pvalues, max_rows = hm_max_rows)
volcanoplot(res.DLNvsHEM.tb, title = "MCMDLN vs. nHEM", max_labels = vp_max_labels,
lfc_limit = vp_lfc_limit)
volcanoplot_interactive(res.DLNvsHEM.tb, title = "MCMDLN vs. nHEM")MCMDLN vs MCM1G
res.DLNvs1G.tb <- res.DLNvs1G %>%
data.frame() %>%
rownames_to_column(var = "ENSEMBL") %>%
as_tibble() %>%
left_join(genedata.filtered, ., by = "ENSEMBL") %>%
mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
res.DLNvs1G.sig <- res.DLNvs1G.tb %>%
filter(threshold == TRUE)
pvalues <- res.DLNvs1G.sig$padj
names(pvalues) <- rownames(res.DLNvs1G.sig)
norm.DLNvs1G.sig <- normalized_counts %>%
data.frame() %>%
dplyr::select(samples$Sample_ID[samples$Class %in% c("MCMDLN", "MCM1G")]) %>%
filter(rownames(.) %in% res.DLNvs1G.sig$ENSEMBL)
write.csv(res.DLNvs1G.tb %>% dplyr::select(-threshold),
file = "results/DESeq/MCMDLN_vs_MCM1G_unfiltered.csv",
row.names = FALSE)
write.csv(res.DLNvs1G.sig %>% dplyr::select(-threshold),
file = "results/DESeq/MCMDLN_vs_MCM1G_filtered.csv",
row.names = FALSE)plot_heatmap(norm.DLNvs1G.sig, hmap_colors, annotation_colors,
genedata.filtered, samples, pvalues, max_rows = hm_max_rows)
volcanoplot(res.DLNvs1G.tb, title = "MCMDLN vs. MCM1G", max_labels = vp_max_labels,
lfc_limit = vp_lfc_limit)